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ABSTRACT 

We explore the evolution of the mass distribution of dust in collision-dominated debris disks, using 
the coUisional code introduced in our previous paper. We analyze the equilibrium distribution and its 
dependence on model parameters by evolving over 100 models to 10 Gyr. With our numerical models, 
we confirm that syste ms reach collision al equilibrium with a mass distribution that is steeper than the 
traditional solution bv lDohnanvil ()1969l ) . Our model yields a quasi steady-state slope of n(m) ^ m^^-*® 



[n{a) 



-3.651 



as a robust solution for a wide range of possible model parameters. We also show that 



a simple power-law function can be an appropriate approximation for the mass distribution of particles 
in certain regimes. The steeper solution has observable effects in the submillimeter and millimeter 
wavelength regimes of the electromagnetic spectrum. We assemble data for nine debris disks that 
have been observed at these wavelengths and, using a simplified absorption efficiency model, show 
that the predicted slope of the particle mass distribution generates SEDs that are in agreement with 
the observed ones. 

Subject headings: methods: numerical - circumstellar matter - planetary systems - infrared: stars 



1. INTRODUCTION 

The total mass within debris disks as well as the in- 
frared excess emission produced by their dust are gen- 
erally calculated assuming the analytic estim ate of the 
distrib ution of masses in the asteroid belt by iDohnan-sdl 
(|1969f ). This solution predicts that the sizes follow a 
power-law, with their numbers increasing with decreas- 
ing size a as n(a) a~^'^ . However, a number of recent 
efforts to model observations of debris disks have found 
it necessary to adopt steeper slopes (jKrist et al.l 120101 : 
IGolimows ki et al.ll2qilD. 

Durda fc Dermottl (l997f l and lO'Brien fc GreenbergI 
120031) have shown numerically and analytically that 
a steep tensile strength curve, i.e., the function 
that gives the minimum energy requ ired to disrupt 

a body catastrophical ly (se e, e.g.. iHolsapple et aTl 

[200l iBenz fc Asphau^ Il999l : iGaspar et al.l re- 
sults in a steeper quasi steady-state distribution than 
the traditional solutio n. Collisional models of the 
dust i n deb ris disks (Thcbault ct al. 2003; Krivov et al.' 
2005t iThebault fc Augereau 20 07; Lohne e t al., ,2008; 
Miiller et aL|[2010t iWvatt et al.H201lh have also shown 



that the dust particles will settle with a size distribu- 
tion slope larger than 3.6, on top of which additional 
structures appear. This steeper distribution has readily 
observable effects at the far-IR and submm wavelengths. 
It also results in higher total dust mass and lower plan- 
etesimal mass estimates for the systems. However, these 
results are often disre garded in observa tional studies. In- 
stead, the traditional iDohnanyil (|1969( ) slope of mass dis- 
tribution is used, likely due to the complexity of the nu- 
merical models, the superposed wavy structures on the 
distributions, and their uncertain regimes of applicabil- 
ity- 

In this paper, we investigate the slope of the mass 
distribution and the physical parameters that influence 
it with the numerical code introduced in IGaspar et al.l 



(2012, hereafter Paper I). Our code calculates the evo- 
lution of the particle mass distribution in collisional sys- 
tems, taking into account both erosive and catastrophic 
collisions. In fj2l we introduce models for the numerical 
analysis of the collisional cascades and give our findings. 
In 331 we analyze the ranges where a simple power-law 
function is an appropriate approximation for the mass 
distribution of particles. In i|4l we compare our model's 
properties to the assumptions of previous analytic solu- 
tions, while in ij5l we generate a set of synthetic spectra 
in order to analyze the effects certain distribution pa- 
rameters have on different parts of the SEDs. In iJSl we 
introduce a simple relation between the Rayleigh-Jeans 
part of the spectral energy distributions and the particle 
size distribution. In fjTl we compare our results to the 
observed far-IR and sub millimeter data for nine sources. 

2. NUMERICAL MODELING 

In this section, we analyze the dust distribution with 
our full numerical code (Paper I). We run a set of nu- 
merical models to study the evolution of the slope of the 
distribution function and its dependence on the model 
parameters. We investigate the time required for the dis- 
tribution to settle into its quasi steady-state and, with a 
wide coverage of the parameter space, we also examine 
the robustness of the solution. 

2.1. Evolution of the reference model 

We set up a reference debris disk as a basis for com- 
parison to all other model runs. This model consists of a 
moderately dense debris disk situated at 25 AU around 
an AO spectral- type star, with a width and height of 
2.5 AU. This radial distance ensures a moderate evolu- 
tion speed. The peak emission is in the mid-infrared, 
with a Rayleigh-Jeans tail in the far-infrared regime, 
which is the primary imaging window for the Herschel 
Space Telescope. The total mass in the debris disk is 
1 , distributed between minimum and maximum par- 
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Table 1 

Numerical, CoUisional, and System parameters of our reference model 



Variable 



Description 



Fiducial value 



System variables 



P 

Mtot 

rio 

R 

AR 

h 

Sp 



Bulk density of particles 

Mass of the smallest particles in the system . . . . 

Mass of the largest particles in the system 

Total mass within the debris ring 

Initial power-law distribution of particle masses 

Distance of the debris ring from the star 

Width of the debris ring 

Height of the debris ring 

Spectral-type of the star 

CoUisional variables 



2.7 g cm ^ 
1.42x10-21 kg 
1.13x10^2 kg 
1 Mg3 

1.87 
25 AU 
2.5 AU 
2.5 AU 
AO 



7 Redistribution power-law 

l3x Power exponent in X particle equation 

a Scaling constant in A/cr 

6 Power-law exponent in Mcr equation 

/m Interpolation boundary for erosive collisions 

fy Fraction of Y/Mcv 

jmax Largest fraction of Y/X at super catastrophic collision boundary 

Qbc Total scaling of the Q* strength curve 

S Scaling of the strength regime of the Q* strength curve 

G Scaling of the gravity regime of the Q* strength curve 

s Power exponent of the strength regime of the Q* strength curve . 

Power exponent of the gravity regime of the Q* strength curve 



TTJW 
1.24 
2.7 X 10-6 
1.23 
10-* 
0.2 
0.5 
1 

3.5 X 10'^ erg/g 
0.3 erg cm^/g2 
-0.38 
1.36 



Numerical parameters 


d 


Neighboring grid point mass ratio 


1.104 





Constant in smoothing weight for large-mass coUisional probability 


lO^mniax 


P 


Exponent in smoothing weight for large-mass coUisional probability 


16 



tide masses that correspond to radii of 5 nm and 1000 
km, when assuming a bulk density of 2.7 g cm"'^. We 
summarize the disk parameters of the reference model in 
Table [H We evolve the reference model for 10 Gyr. 
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Figure 1. Particle mass distribution of the reference model plot- 
ted at various points in time. 



In Figure [U we show the evolution of the particle dis- 
tribution, plotting it at six different points in time up to 
10 Gyr. On the vertical axis we plot logj^p l"^) ^ '^^] ; 
which can be related to the "mass/bin" that is frequently 
used in other simulations. Even though the number den- 
sities decrease with increasing particle masses, the mass 
distribution increases towards the larger masses in this 
representation, as long as the mass distribution slope is 
smaller than 2. 



The smallest particles reach coUisional equilibrium 
first, roughly at 1 Myr, followed by larger particle sizes 
as the system evolves. After 50-100 Myr of evolution, 
the upper, gravity dominated part of the distribution 
(m > 10^'^ kg) also reaches equilibrium. The distribu- 
tion maintains its slope for masses below 10^° kg, which 
roughly corresponds to a planetesimal radius of 100 m. 
The kink in the distribution at the upper end is due to the 
change in the strength cu rve slope (jO'Brien &: Greenbergj 
[200llBottke et al.l]2005l) . 

Structures in the distribution slope, such as waves, 
may in principle occur at the low-mass end when assum- 
ing softer material properties or higher collision veloci- 
ties (see Section [3]). The distribution may also acquire 
some curvature (see Paper I). Because of these effects, 
we evaluate the average slope of the distribution by fit- 
ting a power-law over a large mass range, but one that 
remains below the kink in the distribution. Specifically, 
we fit the slope between masses 10"^ and 10"* kg, which 
roughly correspond to sizes of 0.1 mm and 1 m. 

We next examine the dependence of the quasi steady- 
state distribution slope on the initial conditions , i.e. the 
initial mass-distribution slope rjQ and the initial total 
mass in the disk Mtot- Figure [2] shows the evolution 
of the particle mass distribution and its slope as a func- 
tion of these initial conditions. The left panels show the 
mass distribution after 10 Gyr for different values of the 
two input parameters, while the right panels show the 
evolution of the dust-mass distribution slope. Variations 
in the initial slope, rjo , do not affect the final slope value, 
although the high mass end evolves differently or reaches 
equilibrium at different timescales. A distribution with 
less dust initially (770 < 1-87) also takes more time to 
reach equilibrium. A shallow distribution with an initial 
slope of 770 = 1.5 takes as much as ^ 1 Gyr to reach 
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Figure 2. (Left panels) Particle mass distribution at 10 Gyr, wlien varying the initial mass distribution slope (top) and the total mass of 
the system(bottom). [Right panels) Evolution of the dust-mass distribution slopes when varying the initial mass distribution slope (top) 
and the total mass of the system (bottom). The quasi steady-state distribution slope is practically independent of these initial conditions. 



equilibrium, although such initial d i stribu tion slopes are 
unlikely. As shown bv lLohne et al.l ()2008[ ). the evolution 
of the particle mass distribution is scalable by the total 
mass (or number densities of particles), which is what 
we see in the bottom two panels of Figure [2j All systems 
with different initial masses reach the same equilibrium 
mass distribution, but on different timescales. More mas- 
sive systems evolve on shorter timescales, thus reaching 
their equilibrium more quickly, while less massive sys- 
tems evolve more slowly. 

2.2. The dependence of the quasi steady-state 
distribution function on the collision parameters 

The parameters that describe the outcomes of col- 
lisions, in principle, should be roughly the same for 
all collisional systems. They are the fragmentation 
constants and the par ameters of the strength curve 
()Benz fc Asphaud 119990 . To investigate their effects on 
the evolution of the particle mass distribution, we vary 
their nominal values and evolve the models to the same 
10 Gyr, as we did for the reference model. 

We give here a detailed analysis of the effects of varying 
only five of the twelve parameters (a, b, Qsc, s, S), as 
the remaining seven (7, /3x, /y, /x'"'^ Im, g, G) have 
no significant effects (see Table [T] for the description of 
these parameters). 



In Figure [3l we plot the resulting mass distributions 
when varying the crater ed mass parameters a and h 
(|Koschnv &: Griinll2001airbf ). The parameter a is the to- 
tal scaling and b is the exponent of the projectile's kinetic 
energy in the equation of the cratered mass 

Mer=.(^)V (1) 

As can be seen in Figure 121 the resulting mass distri- 
butions depend on the values of a and b only in the 
gravity dominated regime. At these larger masses, our 
model is incomplete, because we do not include aggre- 
gation. When increasing a, i.e., basically softening the 
materials or increasing the effects of erosions, the num- 
ber of eroded particles in the gravity-dominated regime 
increases rapidly. A similar effect can be observed when 
increasing the value of b. However, within reasonable 
values of a and 6, the variation of the equilibrium par- 
ticle mass distribution slope in the dust mass regime is 
negligible. 

In Figure HI we plot the resulting distributions and 
the evolution of the dust distribution slope when we 
vary the parameters Qsc, S, and s. These are all vari- 
ables in the tensile strength curve, which is given as 
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()Benz fc Asphaudll999l ) 

Q*{a) = lO-^c \S (— )% Gp (—Y] . (2) 

L Vcm/ Vcm/ J erg kg 

The variable Qsc is a global scaling factor, S is the scaling 
of the strength-dominated regime, s is the power depen- 
dence on particle size of the strength-dominated regime, 
G is the scaling of the gravity-dominated regime, and g 
is the power dependence on particle size of the gravity- 
dominated regime. Variations in the gravity-dominated 
regime of the curve (G and g) do not have significant 
effects on the equilibrium dust-mass distribution, so we 
do not consider these parameters further. The tensile 
strength curve has been extensively studied for decades. 
However, as it is dependent on v arious material prop- 
erties and the coUisional velocity (jStewart fc Leinhardtl 
l2009t iLeinhardt fc Stewartil2012t) . its parameters do not 
have universally applicable values. Determining the ten- 
sile strength curve at large and small sizes is also ex- 
tremely difficult experimentally. 

The slope s of the strength curve in the strength- 
dominated regime depends on the Weibull flaw-size dis- 
tribution. Its measured values rang e anywhere between 
-0.7 and -0.3 (jHolsapple et al.|[2002l) . Steeper values of s 
make smaller materials harder to disrupt, which results 
in a steeper dust distribution slope. At s = 1.0, the 
smallest particles are hard enough to resist catastrophic 
disruption even when the projectile mass equals the tar- 
get mass. This results in a mass distribution with a slope 
equal to the redistribution slope 7 = 1.83 at the smallest 
scales, while in the fitted region it is significantly steeper. 
At s = 0.6, the smallest particles are still able to destroy 
each other and generate a dust distribution slope that is 
close to 1.91. 

The scaling constants of the tensile strength curve are 
the dominant parameters in the evolution toward the 
quasi steady-state distribution. When reducing the com- 
plete tensile strength curve scaling Qsc, wave structures 
form more easily, as a particle becomes capable of affect- 
ing the evolution of particles much larger than itself (see 
Paper I). When upscaling the tensile strength curve, the 
quasi steady-state distribution slope starts to resemble 
the redistribution slope, as it is the particle redistribu- 



tions that lead the evolution of the particle mass distri- 
bution. When varying the scaling of only the strength 
side of the curve S, similar effects can be seen. 

We have shown that drastic offsets in the collision pa- 
rameter values result in slight changes to the quasi steady 
state mass distribution slope. We conclude that, for 
all reasonable values of the collisional parameters, the 
quasi steady-state dust-mass distribution slope is larger 
or equal to 1.88. 

2.3. The dependence of the distribution function on 
system variables 

There are a number of parameters that can change 
from one collisional system to another: the material den- 
sity p, the minimum and the maximum particle mass in 
the system mmin and m„iaxj the radial distance i?, height 
h, and width AR of the disk, and the spectral type of 
the central star. All these parameters affect three prop- 
erties of the collisional model: the blow-out mass, the 
collisional velocity, and the number density of particles. 
Varying these parameters will change the timescale of 
the evolution and affect the quasi steady-state distribu- 
tion slope. In this subsection, we analyze the effects of 
varying the radial distance on the equilibrium mass dis- 
tribution through the dependence of the collisional ve- 
locity on radius. Modifying either disk parameters AR 
and h or the spectral-type of the star would have simi- 
lar effects. We defer discussion of the variations in the 
timescales to a future paper. 

In the left panel of Figure [SJ we show the effects of 
varying the radial distance, R, on the mass distribution 
evolved to 10 Gyr. Decreasing the radial distance will 
increase the collisional velocity, resulting in the appear- 
ance of waves at the small-mass end of the distribution. 
It also generates a much more pronounced kink at the 
high mass end. On the other hand, when the velocity 
is decreased at large radii, the low mass end of the dis- 
tribution starts resembling the redistribution function, 
as smaller particles are not destroyed due to the lower 
energy collisions. Moreover, no kink is produced at the 
high mass end. In the right panel of Figure [SI we show 
the evolution of the particle-mass distribution slope as 
a function of the collisional velocity. For high velocity 
collisions, the waves render the fitting of a single mass 
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Figure 4. {Left panels) Particle mass distribution at 10 Gyr, when varying the values of the tensile strength curve parameters Qbc, S and 
s. {Right panels) Evolution of the dust-mass distribution slopes when varying the values of the same parameters. 



distribution slope ill constructed but the underlying slope 
of the wavy mass distribution is slightly steeper than for 
the smaller coUisional velocity case. 

2.4. The dependence of the distribution function on 
numerical parameters 

We discuss here the effects of three non-physical vari- 
ables that appear in the numerical algorithm. They are: 
the mass ratio S between neighboring grid points and the 
parameters of the large particle coUisional cross-section 



smoothing formula, 9 and p (see Paper I). In Figure [SI 
we plot the distribution of the model with varying val- 
ues of (5 at 10 Gyr (left panel) and the evolution in the 
slope of the dust distribution (right panel). The evolu- 
tion of the dust distribution is affected by the number 
of grid points we use, converging at 5 = 1.13; this corre- 
sponds to 801 grid points between our mmin and m^ax 
mass range. Using a lower number of grid points leads 
to errors in the numerical integration for the redistribu- 
tion, leading to an offset larger than 7% for the smallest 
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Figure 5. (Left panel) Particle mass distribution at 10 Gyr, when varying the value of the radial distance of the disk. {Right panel) 
Evolution of the dust-mass distribution slopes when varying the same parameter. 
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Figure 6. (Left panel) Particle mass distribution at 10 Gyr, when varying the resolution of the numerical model. [Right panel) Evolution 
of the dust-mass distribution slopes when varying the same parameter. 



particles in the system, when only using half as many 
grid points. We find that the dust distribution slope is 
practically independent of the smoothing variables Q and 
P- 

2.5. The time to reach quasi steady-state 

In our model calculations, the dust distributions in 
the vast majority of cases reach quasi steady-state by 
10-20 Myr and only in a few cases do they take some- 
what longer. The characteristic time is less than 100 Myr 
for al l realistic cases in agreement with pre vious studies 
(e.g-. lWvatt et al.ll2007l : iLohne et al.ll2008D . This shows 
that, apart from second generation debris disks, the ma- 
jority of debris disks around stars of ages over 100 Myr 
are most likely to be in coUisional quasi steady-state, at 
least for the smallest particles (< 1 cm). However, young 
and e xtended systems, such as P Pic (jSmith fc Terrild 
Il984t iVandenbussche et al.l I2010D . might not be near 
quasi steady-state at the outer parts of their disks, where 
interaction timescales are longer. 

2.6. The robustness of the solution 



One of the most surprising results of the wide range of 
numerical models we computed is the robustness of the 
quasi steady-state distribution. Varying the values of the 
model parameters does not result in significant changes 
in the slope of the distributions. In Figure [71 we show 
the effect on the equilibrium slope of the dust-mass dis- 
tribution function of varying each model parameter. We 
order the parameters on the horizontal axis as a func- 
tion of decreasing magnitude of their effects. Parameters 
that had to be varied by many orders of magnitude to 
have effects on the equilibrium slope are analyzed in log 
space. The plot shows that the dominant parameter, by 
far, is the slope of the strength curve in the strength- 
dominated regime. This is followed by variables that af- 
fect the coUisional velocity {AR and h) and the power b 
of the erosive cratered mass formula (equation [T]). The 
plot also shows that neither of our arbitrarily chosen col- 
lision prescription constants have any significant effects 
on the outcome of the coUisional cascade. The model 
runs predict an equilibrium dust mass distribution slope 
of 7? = 1.88 ± 0.02 (77a = 3.65 ± 0.05), taking the max- 
imum offsets originating from the 50% variation in the 
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Figure 7. Effect on the equilibrium particle mass and size distributions of varying each collisional and system variable by 50%. 



model parameters test as our error. 

3. THE POWER-LAW APPROXIMATION OF THE 
MASS-DISTRIBUTIONS 

Due to the complexity of full numerical models and 
the varying amount of structures in the particle mass- 
distributions presented in the papers detailing them, it 
is still common practice for debris disk SEDs to be cal- 
culated and fitted to observed data assuming a simple 
power-law for the mass-distribution of particles. The 
applicability of a power-law distribution has never been 
probed systematically. Another issue within this sim- 
plification is that the value for the slope of the power- 
law used is generally the tradit i onal 1 1/6 (or 3.5 in size 
space) determined bv lDohnanvil (|1969l ). although this has 
already been proven to be only valid for the constant ma - 
teri al tensile strength case bv iDurda &: DermottI ()1997f ) 
and lO'Brien fc Greenber^ (|2003[). As we have shown in 
the preceding sections, a steeper distribution function 
of 77 = 1.88 ± 0.02 (r/a = 3.65 ± 0.05) is a more likely 
outcome when considering realistic conditions. Similar 
results have been fou nd by IDur da & Dermott ( 1997}; 
lO'Brien fc Greenberal (j2003l ) and IThebault fc Augereaul 
(I2007D . among others. 

In this section, we investigate the extent to which a 
simple power-law is a valid approximation for the mass 
distribution in debris disks. As we are evaluating the 
detectability of effects of deviations from a power-law 
mass distribution on the generated SEDs, we are only 
concerned with particles smaller than ~ 1 mm in size. 
Deviations from a power-law on this scale occur when 
waves appear superposed on the distribution itself, mean- 
ing the conditions that result in the formation of waves 
in the particle mass-distribution need to be studied. 

3.1. Appearance of waves in the particle 
mass- distribution 

Waves appear superposed on the mass-distribution due 
to the discontinuity at the low-mass cut-off. This behav- 
ior has been explained befor e as a result of the build-up 
of particles at the b oundary (|Campo Bagatin et al.lll994 : 
iWvatt et al.ll201ll) . however examining the properties of 
the collisional equation (see Paper I) yields a different 
picture. In Figure IH we plot the d/z integrand of the 



collisional equation 



(X J d^< -TO ^/i ^ (^7713 + ^J.^ J 

+^i-''^Kri + M')^ (5 [X{^l, M) - to] 
+ / dM/z-^'M-" (/i^ +A/3) X 

A{^i,M)H[Y{^JL,M)^m] \ , (3) 



which basically gives the amount of variations produced 
by a /i mass projectile for a certain to mass, either by 
removing to or adding it as an X fragment. In regimes 
marked with a the effective changes in the number 
densities are positive, while in regimes marked with a 
"— " they are negative. With solid black lines we show 
the integrands for the to masses where the first dip of 
the waves are situated for the corresponding velocities. 
It is apparent that for these masses the peaks of the re- 
moval integrands are located exactly at the cut-off. The 
fact that they are peaks is obvious when comparing with 
the integrands of larger masses. If the distributions were 
continuous, this effect would smooth out, but with the 
cut-off located at the peak of the integrand, relatively 
more will be removed of mass to at the first dip than of 
other masses. We also note that the additive peaks of the 
minimum (cut-off) mass, shown with red solid curves, are 
not correlated with the location of the dip in the distri- 
bution (marked with a black solid tick mark) , but shift as 
a function of the collision velocity. At higher collisional 
velocity, where waves are stronger, the largest number of 
cut-off particles are actually produced by the particles in 
the dip, which would provide a negative feedback accord- 
ing to the traditional picture, canceling the waves. Our 
analysis shows that positive feedback is not the dominant 
effect in the production of the waves. 

The location of the peak of the removal term (Term I 
in Paper I) is the determining factor in the formation of 
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Figure 8. The integrands of Equation [3] for selected m masses in system with coUisional velocities of 3 km s~^ (left panel) and 6 km s~^ 
(right panel). The solid black lines give the integrands for the mass at the first wave dip in each system (rriwavo), while the solid red lines 
give the integrands for the minimum mass in the system. The removal peak of the mwavo integrand is at the cut-off, explaining the excess 
removal at the first dip in the waves. The tick marks show the location of the masses in the distribution for which integrands are displayed. 
See text for detailed description of the formation of the wave structure. 



the waves, which can be given by solving 



= a<^ /i-" 



m I /i 3 + J7J, 




M3 



(4) 



where m gives the location of the first dip and M is 
given so that X{fi,M) = m. Unfortunately, an analytic 
solution can not be given, as the derivative is transcen- 
dent, but is easily solvable numerically. The variable 
that determines the solution will be X, which depends 
on the collisional velocity, the tensile strength law, and 
the value of mmin- For a given system, where m^in and 
the tensile strength law do not change as a function of 
radial location, the single variable determining the wavy 
structure is the collisional velocity. The minimum mass 
will change as a function stellar types and some mate- 
rial/velocity dependence of the tensile strength law can 
be expected between systems. 

3.2. Constraining the conditions for the appearance of 
waves in the particle mass-dsitribution 

In Section[2l we have shown that variations in the prox- 
ies for the collisional velocity (such as disk radius) can in 
fact result in wavy size distributions, even within a nar- 
row debris ring. Extended debris disks will be even more 
likely to have higher collisional velocities, as the particles 
with higher /3 values (but less than 0.5) get placed on 
high eccentricity orbits. These small particles from the 
inner rings will collide with the particles in the external 
rings with in creased velocities due to the non-zero colli- 
sional angles (jThebault fc Augereaull2007|) . To study the 
effects of variations only in the collisional velocity, we set 
the collisional velocity directly within our code to specific 
values, not changing any other parameters. We present 
the results from these runs in Figure [9l The figure shows 
that waves start to appear at collisional velocity values 
of 3 km s~^ and above. 

The value of the orbital velocity for a particle on an 
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Figure 9. Particle mass distribution at 10 Gyr, when varying the 
collisional velocity in the system. Waves start to appear at and 
above collisional velocities of 3 kms~^, when constraining all other 
parameters at their fiducial values. Note: The collisional velocity 
was adjusted by hand in these models and not determined by the 
system geometry as described in Paper I, in order to study the 
effects of only variations in the collisional velocity itself. 



elliptical orbit originating from i?in with a certain /3 value 
at i?out radius (where it will collide with a particle on a 
circular orbit) is 



1 - 2/3 



Rn 



(1 - /3) Ru 



(5) 



The angle between the orbital velocities can be expressed 
as 



C — cos 



2i?out (^out — 2a) 



(6) 
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where 



a = Ru 



1 - 2P 



(3 



1-13 

The colhsional velocity is then 



Wcoll 



{Vo 



(vc sin C) 



(7) 
(8) 

(9) 



where Uq is the orbital velocity for the particle on a cir- 
cular orbit at i?out- 

We plot the values of this colhsional velocity for parti- 
cles originating from = 10 AU and 50 AU in systems 
around a solar- and early-type star in Figure [TOl We 
shade the colhsional velocity region above our 3 km s~^ 
limit, where waves start to appear. As we can see, nearby 
rings at AR — 1 AU have low colhsional velocities, re- 
assuring that our prescription for a constant colhsional 
velocity value in narrow debris rings is reasonable. We 
can also see that once we depart from the narrow ring 
assumption {AR > 10 AU), the colhsional velocities do 
increase. 

Around a solar-type star, particles originating from an 
inner radius of 10 AU are able to achieve colhsional ve- 
locities of ~ 3 km s~^; however, the particle size range 
that is able to achieve this limit is very narrow, be- 
tween 0.4 and 0.6 /im. Taking into account the dilu- 
tion of these small particles within a limited narrow size 
range, we conclude that extended debris disks with in- 
ner boundaries outside of 10 - 15 AU around solar type 
stars should be well approximated by a simple power-law 
mass (size) distributio n function . In Paper I, we show a 
compa rison run to Lo hne et al.l (2008) and Wva tt et al.l 
assuming an extended disk between radii of 7.5 
and 15 AU around a solar-type star. The runs by 
iWvatt et all (|20ll and our code sh o w neg ligible / small 
amplitude waves, while iLohne et all (|2008'i show mod- 
erate amplitude waves, likely due to their full three di- 
mensional modeling of the system, taking into account 
smaller particles originating from regions inward of 10 
AU. These results confirm our analytic analysis on the 
conditions required to initiate waves in the dust particle 
mass-distribution. 

We show the same analysis for early-type stars in the 
right panels of Figure [101 The colhsional velocities of 
particles originating from 10 AU with particles on exter- 
nal circular orbits around an AO spectral-type star will 
be significantly larger than our limit for a larger range 
of particle sizes (5 - 10 fim). Particles originating from 
orbits inside of 10 AU will dehnitely initiate waves in ex- 
tended disks or external rings. However, particles origi- 
nating from 50 AU will barely reach colhsional velocities 
at or above our limit and for a limited range of particle 
sizes. We conclude that extended debris disks with inner 
boundaries outside of 50 - 60 AU around early-type stars 
should be well approximated by a simple power-law mass 
(size) distribution function 

Our ID "particle-in-a-box" numerical code adjusts 
the colhsional velocities through the ring thickness and 
height {AR and h), and the resulting orbital inclina- 
tions for the particles. Our results are only moder- 
ately affected over the relatively narrow range of ve- 



locities obtained by varying them (see Figure [T]) . We 
modified the code (see above) to simulate the higher- 
velocity collisions with particles on elliptical orbits to 
determine a threshold for generation of waves in the 
particle distribution. However, our model cannot ex- 
plore the nuances resulting from varying colhsional ve- 
locities in extended systems (such as varying colhsional 
rates and colhsional outcomes), which could also play 
a role in the possible formation of substructures super- 
posed on the generally steeper particle size-distributions 
fsee e.g..lKrivov e t al. 2005; Thcbault & Augcrcau 200'^ 
ILohne et alEooa iMiiller etral.ll2010ll . 

As stated at the end of Section 3.1, our generaliza- 
tion assumed a specific system (with same material prop- 
erties and minimum particle mass). Although we ex- 
plored varying the minimum particle mass when intro- 
ducing different spectral type central stars, both the ten- 
sile strength and the minimum mass may vary with ma- 
terial properties as well. In addition to minerals such as 
basalt, it is likely that grains in the outer zones of debris 
disks contain significant amounts of ice . The mechanical 
properties of icy grains arc explored bv lBenz fc Asphaugl 
(1999) and Leinhardt & Stewart (200i), yielding roughly 
a factor of 3 in difference between the tensile strength of 
basalt and ice. Given the relative insensitivity of our 
results to the scaling of the grain strength curve {Qsc - 
see Figure [7]) , these properties are similar enough that 
our basic conclusions about the grain size distribution 
and region of applicability of the power law approxima- 
tion to it are essentially unchanged even outside of the 
iceline. 

4. COMPARISON TO PREVIOUS ANALYTIC APPROACHES 

The slope of the size distribution slope in colhsional 
cascades has been investigated with analytic approaches 
wi th increasi n g com plexity since the pioneering work 
of iDohnanvil ()1969D . His assumptions of self-similar 
fragmentation, constant material strength, and collision 
rates proportional to the geometric cross section lead to 
the well known and now widely used steady state size- 
distribution slope of rjm = 11/6 {rja = 3.5). The validity 
of his assumption s has been inyestigat ed in a number of 
analytic studies. iTanaka et al.l ()1996[ ) found that vary- 
ing the power of the colhsional cross section {m'^ , where 
u = 2/3 for a simple geometric cross section) will set 
the steady state distribution slope as rjm = {i^ + 3)/2. 
Such non-geometric cross sections (collision rates) can 
be the results of mass weighted collision probabilities, 
varying amounts of g rain porosity, mass dependen t colh- 
sional velocities, etc. lO'Brien fc Greenbergj (|2003D inves- 
tigated whether a non-constant tensile strength value, 
defined by a power-law with slope of s, will affect the 
steady state particle size distribution. They found a lin- 
ear relationship between the two power-law slopes, yield- 
ing a steady state solut ion produc i ng m ore small parti- 
cles than the traditional IDohnanvil (|1969f ) solution, when 
considering realistic tensile stren gth slopes. An ana- 
lytic solution similar to t he one bv O'Brie n fc Greenbergl 
(2003) was presented in IWvatt e t al. ( 20'TTI ). who also 
analyzed their equations numerically, revealing sub- 
structures superimposed on t he po wer-law size distri- 
butions. iBelvaev fc Rafiko^ ()201lD expanded on the 
lO'Brien fc Greenbera ()2003i) model by allowing a vari- 
able largest fragment mass to target mass ratio {X/M 
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Figure 10. The coUisional velocity between particles originating from a radius of iJin and colliding with particles at radius -Rout, with 
eccentricities calculated from their /3 values. We show cases for a solar-type star and an AO spectral-type star, with values of 10 and 
50 AU. 



Table 2 

Comparing the properties of analytic solutions with our numerical model 



Model 



Collision rate 



XjM 



Redistribution 



d In Ti(m) 
rf In m 



Dohnanvi f 19691 
Tanaka et al. ri996') 
O'Brien fc Greenberg (2003) 
Bel yaev fc Rafikov (20111 

This work 



2/3 



my 



Constant 
Constant 
Constant 
Variable 



Constant 
Constant 
power-law 
power-law 



power-law 
power-law 
power-law 
power-law 



variable 



Variable 



power-law power-law 



11/6 

{ly + 3) /2 

(lH-s)/(6 + s) 

slowly varying function of m 

I'unction of model parameters, 

largely robust around -1.88 



in our model), follo wing the iFuiiwara et alJ (1197711 ex- 
perimental results. iBelvaev fc Rafikovl ()2011l1 found a 
particle distribution that could not be described by a 
single power-law, rather by a slowly varying, mass de- 
pendent power-law function. We summarize the charac- 
teristic properties of these analytic models in Table [21 
A common issue with all of these analytic solutions is 
that they assume a steady state system, meaning that 
the mass flux through mass m is zero [dF{m)/dm = 0]. 
In reality, without mass input at the high mass end the 
systems lose mass and decay, meaning the mass flux is 
never zero in a system, but rather is a time and mass 
dependent variable. The solution of the cascade is ac- 
tually an eigenvalue problem, which is difficult to solve 



analytically. 

Our numerical code allows investigating the quasi 
steady state [dF{m)/dm ^ 0] solution of coUisional 
cascades, without the aforeme ntioned simplifica .tions. 
In Figure [TTl we c ompare the IDohnanvH ([l96^ and 
lO'Brien fc Greenberg (2003 ) steady state mass distribu- 
tion slope results with the ones given by our code at 
varying values of tensile stren gth slope. At s = (con- 
stant Q*j2,) the lO'Brien fc Gr eenberg (20^1) solution is 
equal to the traditional iDohnanvil ^1969f pvalue. but it 
gives steeper mass-distribution slopes for lower values of 
s. Similar behavior is seen with our code; however, at 
s = our prediction for the quasi steady state slope is 
rjm = 1.85 {ria — 3.55), which is still steeper than the 



Modeling CoUisional Cascades: Steep Dust-Size Distributions 



11 



E 




-1 -0.8 -0.6 -0.4 -0.2 
s 



Figure 11. Comparing the particle rnass-distribution 

slope vs. tensile-s trength curve slope of Doh nanvii tl969i ). 
lO'Brien fc Greenberei Il2003i) . and our current work. 

traditional value. At much lower values of s (< —0.25) 
we predict a mass distrib ution slope that is less steep 
than the ones given by the lO'Brien fc Greenber^ (|2003f ) 
formula, while our results intersect at s = —0.25, yield- 
ing a quasi steady state slope of rjm — 1.87 {rja = 3.61). 
Since the point of this intersection is close to the theoret- 
ical and experimental value of the tensile strength curve 
slope, the predictions from a pure steady state model and 
a quasi steady state model will agree within errors. This 
renders the testing of different models difficult. How- 
ever, our numerical calculations provide a general new 
estimate for the particle size distribution slopes with sig- 
nificantly more complete physics compared to previous 
analytic approaches. 

5. SYNTHETIC SPECTRA 

In the following sections, we compare the emission that 
results from the predicted particle-mass distributions to 
observations. As a first step, we generate an array of syn- 
thetic spectra using realistic astronomical silicate emis- 
sion properties. We then analyze how the spectra are 
influenced by the particle mass distribution function. 

The flux emitted by a distribution of particle masses 
at a certain frequency is equal to 

tlTT 

^- = J^ J da n(a)a'Qabs (a, i^) (T) , (10) 

where Qabs is the absorption efficiency coefficient, Bi, (T) 
is the blackbody function, and is the total volume of 
the emitting region. Since in infrared astronomy it is 
customary to express the flux density as a function of 
wavelength, we rewrite this also as 

= ^ / d« n(a)a2gabs (a, A) Bx (T) . (11) 

The exact function of the absorption efficiencies of par- 
ticles in the interstellar medium or in circumstellar disks 
is largely unknown. The most commonly used particle 
types for SED calculations are artificial astronomical sil- 
icate (the properties of which are adjusted to reproduce 
the typical 10 fim silicate feature and measured labo- 
rator y dielectric functions) and graphite ()Draine fc Led 
11984) . In Figure [lU we plot the absorption efficiency 
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Figure 12. Absorption efficiency of astronomical silicates as a 
function of wavelength for a range of particle sizes between 0.1 and 
80 /xm. Solid lines are the ~ A"'' approximations to the long- 
wavelength regimes of the curves that we employ in this work. 

as a function of wavelength for a few astronomical sil- 
icate particle sizes. Particles larger than 10 fj,m have 
nearly constant absorption efficiency curves at shorter 
wavelengths (A < 27ra, where a is the particle radius) 
with Qabs = 1, which is followed by a power- law cut off. 
The slope of this power-law becomes constant for wave- 
lengths larger than ^ Stto, and is commonly denoted by 
the variable /?. Astron omical silicates of a ll sizes have a 
typical value of /? = 2 (|Draine fc Led[l98l . 

In Figure [131 we show synthetic SEDs, all scaled to the 
same flux level at 1000 fim. The top panel shows spectra 
calculated around an AO spectral-type star, with debris 
rings placed at various distances between 2 and 292 AU. 
The minimum particle size cut-off was set at ~ 5 /um, in 
accordance with our model (see Paper I). All disks with 
radial distances below ^ 130 AU have a common slope 
for wavelengths larger than 250 /xm, and the furthest disk 
at 292 AU joins this common slope around ~ 350 /im. 

The blow-out size in a system depends on grain struc- 
ture (porosity) and the exact value of the optical con- 
stants for small grains (which is largely unknown and is 
a function of grain material). For this reason, we also 
calculated synthetic SEDs for a debris ring at 25 AU 
around an AO spectral-type star, with the minimum par- 
ticle size of the distribution artificially cut off at sizes 
between 5 and 30 /im. (Note that we normally calculate 
the blow-out mass self-consistently as described in Pa- 
per I.) We plot these SEDs in the middle panel of Figure 
[T3l The offsets between the SEDs become apparent for 
wavelengths shorter than 200 /im, while for longer wave- 
lengths, the emission profiles agree and have a common 
pseudo Rayleigh- Jeans slope. 

Finally, we explore the dependence of the SED on the 
slope of the quasi steady-state particle-mass distribu- 
tion. The bottom panel of Figure [13] shows synthetic 
SEDs generated for a debris ring at 25 AU around an 
AO spectral type star, with a minimum particle cut-off 
size at 5 /im, but with particle mass distribution slopes 
between 1.81 and 1.99. These plots show that the slope 
of the Rayleigh- Jeans part of the emission is greatly in- 
fluenced by the particle size distribution slope. In fact, 
it depends almost solely on this slope, with the temper- 
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Figure 13. Synthetic SEDs for an array of model power-law par- 
ticle mass distributions, with varied parameters. The fluxes are 
scaled to match at 1000 /jm. In the top panel, we show the syn- 
thetic SEDs generated for a variety of radial distances; in the mid- 
dle panel, we show synthetic SEDs generated for a variety of min- 
imum cut-offs (1, 2, 4, 8, 15 and 30 /^m); and in the bottom panel, 
we show synthetic SEDs generated for a variety of particle-mass 
distribution slopes (1.80, 1.84, 1.88, 1.92, 1.96 and 1.99). Note: In 
the middle-panel we vary the minimum cut-off of the particle mass 
distributions, even though it is a parameter inherently set by equa- 
tions in our coUisional model. Since in reality, the placement of the 
cut-off is set by the optical properties and structural build of the 
micron size particles, it is generally treated as a variable in SED 
models. In the plot we show that such variations in the placement 
of the minimum cut-off do not affect the long-wavelength part of 
the SEDs. 



ature of the grains having mild effects at large orbital 



distances. We al so performed our tes ts with dirty-ice 
optical constants ([Preibisch et al.|[T993D and found very 
similar results, showing that our results are also indepen- 
dent on particle types assumed. 

6. RELATION BETWEEN THE PARTICLE MASS 
DISTRIBUTION AND THE SED 

The absorption efficiency curves can be simplified and 
described as 

[ 1 A < 47ra 

Qabs(A,a)oc|^^^, A>4™ 

where x is a scaling constant for the power-law part of 
the function. Fitting the silicate absorption efficiency 
functions, we find 



12 



-0.5 



(12) 



Using this simplified absorption efficiency model and as- 
suming that all particles contribute to the Rayleigh- Jeans 
tail of the SED with their own Rayleigh- Jeans emission, 
we estimate the emitted flux density at long wavelengths 
as 



2p7rkbTCdisk 
D2A2 



da a^-"" 



10" 



da 



(13) 



Here we assumed a /? parameter that is independent of 
the particle size. The variable Cdisk is the number density 
scaling (see Paper I), kb is the Boltzmann constant, T is 
the temperature of the dust grains (which we also assume 
to be particle size independent), and £) is the distance 
of the system from the observer. The quantity 77^ is the 
quasi steady-state particle size distribution slope, and 
can be calculated from the mass distribution slope as 
rja — Srj — 2. Integrating these functions, we get 



F,(A) = Ci A^- 



+ C2 A 



(14) 



where 



Ci = X ^-^^ — (15) 



C2 = 



2t,^kbTCdi.k (4^)"""' 



(16) 



Assuming /3 — 2, which is appropriate for astronomical 
silicates, we find that the slope of the SED is equal to —rja 
for the short wavelength part of the Rayleigh- Jeans tail 
of the SED and 1 — 77a for the long wavelength regime, 
where (5abs(A,a) = 1 for the particles contributing the 
most to the em i ssion. Similar results have been found by 
iWvatt k Den^ ()2002 ^. This behavior explains why the 
slope of the flux density in the long wavelength regime 
is not dependent on the optical properties of the grains, 
as Qabs(A, a) = 1 for all grain types at these wavelengths 
that effectively contribute to it. 

Our models yield a quasi steady-state distribution 
slope of rja ~ 3.65, meaning that the Rayleigh- Jeans tail 
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end of the SEDs should be proportional to 

OC A-2-65 ^ (17) 

as long as the particles are in collisional quasi steady- 
state. 

7. COMPARISON TO OBSERVATIONS 

To compare the computed spectra of quasi steady- 
state collisional disks to data, we assembled the available 
data for debris disks with far-IR and submillimeter ob- 
servations. As a result of our analysis in ^ where we 
determined the wavelength range that is least sensitive 
to parameters, we use only data at wavelengths larger 
than 250 /xm. To fit a power-law to the Rayleigh- Jeans 
regime of the SEDs, we need a minimum of three data 
points above our wavelength cut-off. We found a total 
of only nine sources that fulfill these requirements. We 
present the far-IR/submillimeter fluxes for these sources 
in Table [3l Occasionally, published submillimeter mea- 
surements do not account for systematic errors. In these 
cases, we applied a total of 30% error to all ground based 
measurements at 350 and 450 /im and 15% for all Her- 
schel data and measurements above 850 /im. We also 
made sure that the data included all the flux from each 
source and applied an aperture correction estimate oth- 
erwise. All corrections are listed as notes in Table [31 

We perform individual power-law fits to the data of 
each source as well as a fit to all sources simultaneously 
with a common Rayleigh- Jeans slope. In Figure I14[ we 
present the photosphere subtracted excess emissions for 
each source in the left panels and plot the best-fit power- 
law spectrum of the form 

F^:^Ax( ^ , (18) 

V200 ^my ' ^ ' 

obtained from individual fits. We show in the right pan- 
els of Figure [14] the error contours of the slope and nor- 
malization of the power-law at the 1, 2, and 3ct levels. 
The plots also indicate the Xmin/d-o.f. (the minimum of 
the reduced x^) of each fit. The solid red line represents 
the R ayleigh- Jeans slope calculated from the iDohnanyil 
(|1969f ) analytic solution, the green band represents the 
best slope given by our reference model calculation (in- 
cluding errors from 50% variations in the slope of the 
strength curve, see Figure [7|), and the blue band yields 
our global fit solution of 

/ = 2.58 ±0.06. (19) 

The global fit and our reference model agree within the 
errors of the prediction. 

8. CONCLUSIONS 

In this paper, we used our numerical model introduced 
in Paper I to follow the evolution of a distribution of 
particle masses. Our numerical model has been built to 
ensure mass conservation and that the resulting distribu- 
tion of particles is not artificially offset due to numerical 
errors, as the integrations of the model span over 40 or- 
ders of magnitude in mass. In §2.41 of this paper, we 
demonstrate that lower precision integrations can lead 
to shallower particle distributions. 

We varied all twelve collisional, all six system, and all 
three numerical variables of our model and examined the 



effects of these variations on the evolution of the parti- 
cle mass distribution. The quasi steady-state particle 
distribution of the collisional system is extremely robust 
against variations in its variables, with the strongest ef- 
fects occurring from ch a nges to th e tensile strength curve 
tHolsapple_et_alJ 120021 : iBenz fc Asphaug 1999). Even 
these variations have mild effects on the slope of the par- 
ticle mass distribution, modifying it only between the 
values of 1.84 and 1.94 (3.52 and 3.82 in size space, re- 
spectively). We find the dust mass distribution of our 
reference model to be 1.88 (3.65 in size space). We find 
that waves occur when the collisional velocities are high 
or when particle strengths are low at the mass distribu- 
tion cut-off, where the radiation force blowout dominates 
the dynamics. Nonetheless, in ij3l we show that a sim- 
ple power-law mass (size) distribution is an appropriate 
approximation even for extended disks with components 
only outside of 10 - 15 AU for solar- and 50 - 60 AU for 
early- type stars. 

The Rayleigh-Jeans tail of the debris disk SEDs is dom- 
inated by the medium sized particles, whose mass distri- 
bution is less affected by possible wavy structures. We 
derive a simple formula that gives the slope of the mea- 
sured flux density in the Rayleigh-Jeans part of the SEDs 
as 

This implies that the mass distribution slope itself could, 
in principle, be measured from long-wavelength obser- 
vations. We assemble a list of nine debris disks that 
have been measured at the far-IR, submillimeter, and 
millimeter wavelengths and examine the Rayleigh-Jeans 
slope of their emissions. Our predictions of a slope of 
/ = 2.65 ± 0.05 agrees well with the observations, which 
have a global slope flt of / = 2.58 ± 0.06. 

Support for this work was provided by NASA through 
Contract Number 1255094 issued by JPL/Cahech. We 
thank Viktor Zubko and Karl Misselt for providing the 
numerical code to calculate the optical properties of large 
grains. 
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Observational data of debris disks 
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Figure 14. Observed SEDs of debris disks with submillimeter and millimeter data. The left panels are the photosphere-subtracted fluxes 
of the excess emissions with the best fitting slopes, while t he right pa nels are the 68%, 95% and 99% confidence contours of the individual 
fits. The error contours also show the slope given by the IDohnanvil 11969) mass distribution function in red, the value predicted by our 
numerical code in the green band, and the best global fit with errors in the blue band. 



